{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_helper": "subslide_end",
     "slide_type": "subslide"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Reading netCDF data\n",
    "- requires [numpy](http://numpy.scipy.org) and netCDF/HDF5 C libraries.\n",
    "- Github site: https://github.com/Unidata/netcdf4-python\n",
    "- Online docs: http://unidata.github.io/netcdf4-python/\n",
    "- Based on Konrad Hinsen's old [Scientific.IO.NetCDF](http://dirac.cnrs-orleans.fr/plone/software/scientificpython/) API, with lots of added netcdf version 4 features.\n",
    "- Developed by Jeff Whitaker at NOAA, with many contributions from users."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Interactively exploring a netCDF File\n",
    "\n",
    "Let's explore a netCDF file from the *Atlantic Real-Time Ocean Forecast System*\n",
    "\n",
    "first, import netcdf4-python and numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "internals": {
     "frag_number": 2,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import netCDF4\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 2,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Create a netCDF4.Dataset object\n",
    "- **`f`** is a `Dataset` object, representing an open netCDF file.\n",
    "- printing the object gives you summary information, similar to *`ncdump -h`*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 4,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'netCDF4._netCDF4.Dataset'>\n",
      "root group (NETCDF4_CLASSIC data model, file format HDF5):\n",
      "    Conventions: CF-1.0\n",
      "    title: HYCOM ATLb2.00\n",
      "    institution: National Centers for Environmental Prediction\n",
      "    source: HYCOM archive file\n",
      "    experiment: 90.9\n",
      "    history: archv2ncdf3z\n",
      "    dimensions(sizes): MT(1), Y(850), X(712), Depth(10)\n",
      "    variables(dimensions): float64 MT(MT), float64 Date(MT), float32 Depth(Depth), int32 Y(Y), int32 X(X), float32 Latitude(Y, X), float32 Longitude(Y, X), float32 u(MT, Depth, Y, X), float32 v(MT, Depth, Y, X), float32 temperature(MT, Depth, Y, X), float32 salinity(MT, Depth, Y, X)\n",
      "    groups: \n"
     ]
    }
   ],
   "source": [
    "f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc')\n",
    "print(f) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 4,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Access a netCDF variable\n",
    "- variable objects stored by name in **`variables`** dict.\n",
    "- print the variable yields summary info (including all the attributes).\n",
    "- no actual data read yet (just have a reference to the variable object with metadata)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 6,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['MT', 'Date', 'Depth', 'Y', 'X', 'Latitude', 'Longitude', 'u', 'v', 'temperature', 'salinity'])\n",
      "<class 'netCDF4._netCDF4.Variable'>\n",
      "float32 temperature(MT, Depth, Y, X)\n",
      "    coordinates: Longitude Latitude Date\n",
      "    standard_name: sea_water_potential_temperature\n",
      "    units: degC\n",
      "    _FillValue: 1.2676506e+30\n",
      "    valid_range: [-5.078603  11.1498995]\n",
      "    long_name:   temp [90.9H]\n",
      "unlimited dimensions: MT\n",
      "current shape = (1, 10, 850, 712)\n",
      "filling on\n"
     ]
    }
   ],
   "source": [
    "print(f.variables.keys()) # get all variable names\n",
    "temp = f.variables['temperature']  # temperature variable\n",
    "print(temp) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 6,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## List the Dimensions\n",
    "\n",
    "- All variables in a netCDF file have an associated shape, specified by a list of dimensions.\n",
    "- Let's list all the dimensions in this netCDF file.\n",
    "- Note that the **`MT`** dimension is special (*`unlimited`*), which means it can be appended to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 8
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('MT', <class 'netCDF4._netCDF4.Dimension'> (unlimited): name = 'MT', size = 1)\n",
      "('Y', <class 'netCDF4._netCDF4.Dimension'>: name = 'Y', size = 850)\n",
      "('X', <class 'netCDF4._netCDF4.Dimension'>: name = 'X', size = 712)\n",
      "('Depth', <class 'netCDF4._netCDF4.Dimension'>: name = 'Depth', size = 10)\n"
     ]
    }
   ],
   "source": [
    "for d in f.dimensions.items():\n",
    "    print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 9
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Each variable has a **`dimensions`** and a **`shape`** attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 10
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('MT', 'Depth', 'Y', 'X')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temp.dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 11,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 10, 850, 712)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temp.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 11,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Each dimension typically has a variable associated with it (called a *coordinate* variable).\n",
    "- *Coordinate variables* are 1D variables that have the same name as dimensions.\n",
    "- Coordinate variables and *auxiliary coordinate variables* (named by the *coordinates* attribute) locate values in time and space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 13,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'netCDF4._netCDF4.Variable'>\n",
      "float64 MT(MT)\n",
      "    long_name: time\n",
      "    units: days since 1900-12-31 00:00:00\n",
      "    calendar: standard\n",
      "    axis: T\n",
      "unlimited dimensions: MT\n",
      "current shape = (1,)\n",
      "filling on, default _FillValue of 9.969209968386869e+36 used\n",
      "<class 'netCDF4._netCDF4.Variable'>\n",
      "int32 X(X)\n",
      "    point_spacing: even\n",
      "    axis: X\n",
      "unlimited dimensions: \n",
      "current shape = (712,)\n",
      "filling on, default _FillValue of -2147483647 used\n"
     ]
    }
   ],
   "source": [
    "mt = f.variables['MT']\n",
    "depth = f.variables['Depth']\n",
    "x,y = f.variables['X'], f.variables['Y']\n",
    "print(mt)\n",
    "print(x)                 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 13,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Accessing data from a netCDF variable object\n",
    "\n",
    "- netCDF variables objects behave much like numpy arrays.\n",
    "- slicing a netCDF variable object returns a numpy array with the data.\n",
    "- Boolean array and integer sequence indexing behaves differently for netCDF variables than for numpy arrays. Only 1-d boolean arrays and integer sequences are allowed, and these indices work independently along each dimension (similar to the way vector subscripts work in fortran)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 15
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[41023.25]\n"
     ]
    }
   ],
   "source": [
    "time = mt[:]  # Reads the netCDF variable MT, array of one element\n",
    "print(time) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 16
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[   0.  100.  200.  400.  700. 1000. 2000. 3000. 4000. 5000.]\n"
     ]
    }
   ],
   "source": [
    "dpth = depth[:] # examine depth array\n",
    "print(dpth) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 17,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape of temp variable: (1, 10, 850, 712)\n",
      "shape of temp slice: (6, 425, 356)\n"
     ]
    }
   ],
   "source": [
    "xx,yy = x[:],y[:]\n",
    "print('shape of temp variable: %s' % repr(temp.shape))\n",
    "tempslice = temp[0, dpth > 400, yy > yy.max()/2, xx > xx.max()/2]\n",
    "print('shape of temp slice: %s' % repr(tempslice.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 17,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What is the sea surface temperature and salinity at 50N, 140W?\n",
    "### Finding the latitude and longitude indices of 50N, 140W\n",
    "\n",
    "- The `X` and `Y` dimensions don't look like longitudes and latitudes\n",
    "- Use the auxilary coordinate variables named in the `coordinates` variable attribute, `Latitude` and `Longitude`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 19
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'netCDF4._netCDF4.Variable'>\n",
      "float32 Latitude(Y, X)\n",
      "    standard_name: latitude\n",
      "    units: degrees_north\n",
      "unlimited dimensions: \n",
      "current shape = (850, 712)\n",
      "filling on, default _FillValue of 9.969209968386869e+36 used\n"
     ]
    }
   ],
   "source": [
    "lat, lon = f.variables['Latitude'], f.variables['Longitude']\n",
    "print(lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 20,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Aha!  So we need to find array indices `iy` and `ix` such that `Latitude[iy, ix]` is close to 50.0 and `Longitude[iy, ix]` is close to -140.0 ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 20,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# extract lat/lon values (in degrees) to numpy arrays\n",
    "latvals = lat[:]; lonvals = lon[:] \n",
    "# a function to find the index of the point closest pt\n",
    "# (in squared distance) to give lat/lon value.\n",
    "def getclosest_ij(lats,lons,latpt,lonpt):\n",
    "    # find squared distance of every point on grid\n",
    "    dist_sq = (lats-latpt)**2 + (lons-lonpt)**2  \n",
    "    # 1D index of minimum dist_sq element\n",
    "    minindex_flattened = dist_sq.argmin()    \n",
    "    # Get 2D index for latvals and lonvals arrays from 1D index\n",
    "    return np.unravel_index(minindex_flattened, lats.shape)\n",
    "iy_min, ix_min = getclosest_ij(latvals, lonvals, 50., -140)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 22
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Now we have all the information we need to find our answer.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 23
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "```\n",
    "|----------+--------|\n",
    "| Variable |  Index |\n",
    "|----------+--------|\n",
    "| MT       |      0 |\n",
    "| Depth    |      0 |\n",
    "| Y        | iy_min |\n",
    "| X        | ix_min |\n",
    "|----------+--------|\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 24
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### What is the sea surface temperature and salinity at the specified point?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 25,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 6.4631 degC\n",
      "32.6572 psu\n"
     ]
    }
   ],
   "source": [
    "sal = f.variables['salinity']\n",
    "# Read values out of the netCDF file for temperature and salinity\n",
    "print('%7.4f %s' % (temp[0,0,iy_min,ix_min], temp.units))\n",
    "print('%7.4f %s' % (sal[0,0,iy_min,ix_min], sal.units))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 25,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Remote data access via openDAP\n",
    "\n",
    "- Remote data can be accessed seamlessly with the netcdf4-python API\n",
    "- Access happens via the DAP protocol and DAP servers, such as TDS.\n",
    "- many formats supported, like GRIB, are supported \"under the hood\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 27
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The following example showcases some nice netCDF features:\n",
    "\n",
    "1. We are seamlessly accessing **remote** data, from a TDS server.\n",
    "2. We are seamlessly accessing **GRIB2** data, as if it were netCDF data.\n",
    "3. We are generating **metadata** on-the-fly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 28,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p5deg/GFS_Global_0p5deg_20230525_1200.grib2/GC\n"
     ]
    }
   ],
   "source": [
    "import datetime\n",
    "date = datetime.datetime.now()\n",
    "# build URL for latest synoptic analysis time\n",
    "URL = 'https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p5deg/GFS_Global_0p5deg_%04i%02i%02i_%02i%02i.grib2/GC' %\\\n",
    "(date.year,date.month,date.day,6*(date.hour//6),0)\n",
    "# keep moving back 6 hours until a valid URL found\n",
    "validURL = False; ncount = 0\n",
    "while (not validURL and ncount < 10):\n",
    "    print(URL)\n",
    "    try:\n",
    "        gfs = netCDF4.Dataset(URL)\n",
    "        validURL = True\n",
    "    except RuntimeError:\n",
    "        date -= datetime.timedelta(hours=6)\n",
    "        ncount += 1       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 28,
     "slide_helper": "subslide_end",
     "slide_type": "subslide"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'netCDF4._netCDF4.Variable'>\n",
      "float32 Temperature_surface(time1, lat, lon)\n",
      "    long_name: Temperature @ Ground or water surface\n",
      "    units: K\n",
      "    abbreviation: TMP\n",
      "    missing_value: nan\n",
      "    grid_mapping: LatLon_Projection\n",
      "    coordinates: reftime time1 lat lon \n",
      "    Grib_Variable_Id: VAR_0-0-0_L1\n",
      "    Grib2_Parameter: [0 0 0]\n",
      "    Grib2_Parameter_Discipline: Meteorological products\n",
      "    Grib2_Parameter_Category: Temperature\n",
      "    Grib2_Parameter_Name: Temperature\n",
      "    Grib2_Level_Type: 1\n",
      "    Grib2_Level_Desc: Ground or water surface\n",
      "    Grib2_Generating_Process_Type: Forecast\n",
      "    Grib2_Statistical_Process_Type: UnknownStatType--1\n",
      "unlimited dimensions: \n",
      "current shape = (129, 361, 720)\n",
      "filling off\n",
      "<class 'netCDF4._netCDF4.Variable'>\n",
      "float64 time1(time1)\n",
      "    units: Hour since 2023-05-25T12:00:00Z\n",
      "    standard_name: time\n",
      "    long_name: GRIB forecast or observation time\n",
      "    calendar: proleptic_gregorian\n",
      "    _CoordinateAxisType: Time\n",
      "unlimited dimensions: \n",
      "current shape = (129,)\n",
      "filling off\n",
      "<class 'netCDF4._netCDF4.Variable'>\n",
      "float32 lat(lat)\n",
      "    units: degrees_north\n",
      "    _CoordinateAxisType: Lat\n",
      "unlimited dimensions: \n",
      "current shape = (361,)\n",
      "filling off\n",
      "<class 'netCDF4._netCDF4.Variable'>\n",
      "float32 lon(lon)\n",
      "    units: degrees_east\n",
      "    _CoordinateAxisType: Lon\n",
      "unlimited dimensions: \n",
      "current shape = (720,)\n",
      "filling off\n"
     ]
    }
   ],
   "source": [
    "# Look at metadata for a specific variable\n",
    "# gfs.variables.keys() will show all available variables.\n",
    "sfctmp = gfs.variables['Temperature_surface']\n",
    "# get info about sfctmp\n",
    "print(sfctmp)\n",
    "# print coord vars associated with this variable\n",
    "for dname in sfctmp.dimensions:   \n",
    "    print(gfs.variables[dname])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 28,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##Missing values\n",
    "- when `data == var.missing_value` somewhere, a masked array is returned.\n",
    "- illustrate with soil moisture data (only defined over land)\n",
    "- white areas on plot are masked values over water."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 31
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape=(361, 720), type=<class 'numpy.ma.core.MaskedArray'>, missing_value=nan\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmwElEQVR4nO29f3RV1Z33/45AIkLIN0jJDwlMpJhWA5kOqECdGhWQtOIoXaOdzlCtTpe2CERgyaDrabGrA+KMRKgtz7SPBcE6+MwSHJ0CEsYGx4fHpxo1As7KWMXyQ9KsakjAoYnA+f4R9mXfk/Njn3P2+XHveb/Wugtyz7nnnLvv/vHen/35fHaBYRgGCCGEEEISxAVxPwAhhBBCiBkKFEIIIYQkDgoUQgghhCQOChRCCCGEJA4KFEIIIYQkDgoUQgghhCQOChRCCCGEJA4KFEIIIYQkjsFxP4Afzp49i48++gjFxcUoKCiI+3EIIYQQooBhGDhx4gQqKytxwQXONpKcFCgfffQRqqqq4n4MQgghhPjg8OHDGDNmjOM5OSlQiouLAfR/wREjRsT8NIQQQghRoaenB1VVVZlx3ImcFChiWWfEiBEUKIQQQkiOoeKeQSdZQgghhCQOChRCCCGEJA4KFEIIIYQkDk8CZf369Zg0aVLG92PatGnYsWNH5vidd96JgoKCrNfUqVOzrtHb24sFCxZg1KhRGDZsGG6++WYcOXJEz7chhBBCSF7gSaCMGTMGjzzyCN544w288cYbuP766/EXf/EXOHDgQOac2bNn49ixY5nX9u3bs67R2NiIbdu2YcuWLXj11Vdx8uRJ3HTTTThz5oyeb0QIIYSQnKfAMAwjyAVGjhyJf/iHf8Ddd9+NO++8E8ePH8fzzz9veW53dzc+97nPYfPmzbj99tsBnM9psn37dtx4441K9+zp6UFJSQm6u7sZxUMIIYTkCF7Gb98+KGfOnMGWLVvw6aefYtq0aZn3W1paMHr0aFx22WX4zne+g87Ozsyx1tZWfPbZZ5g1a1bmvcrKStTW1mLv3r229+rt7UVPT0/WixBCCCH5i2eBsm/fPgwfPhxFRUW49957sW3bNlx++eUAgIaGBvzyl7/Eyy+/jMceewyvv/46rr/+evT29gIAOjo6UFhYiNLS0qxrlpWVoaOjw/aeq1atQklJSebFLLKEEEJIfuM5UVtNTQ3efvttHD9+HM899xzuuOMO7NmzB5dffnlm2QYAamtrMWXKFIwbNw6/+tWvMHfuXNtrGobhmLRl+fLlWLx4ceZvkYmOEEIIIfmJZ4FSWFiIz3/+8wCAKVOm4PXXX8fatWvxT//0TwPOraiowLhx4/Dee+8BAMrLy9HX14eurq4sK0pnZyemT59ue8+ioiIUFRV5fVRCCCGE5CiBU90bhpFZwjHz8ccf4/Dhw6ioqAAATJ48GUOGDEFzczNuu+02AMCxY8ewf/9+PProo0EfhRBCCMlb2g5lrxzUjT1s+b44Jt4X5+UangTKgw8+iIaGBlRVVeHEiRPYsmULWlpasHPnTpw8eRIrVqzA17/+dVRUVODDDz/Egw8+iFGjRuHWW28FAJSUlODuu+/GkiVLcPHFF2PkyJFYunQpJk6ciBkzZoTyBQkhhKQLeWC2Grw3d03H1j1XYe61v8G80r2xDuBWzwc4iw+v181VoeJJoPz+97/HvHnzcOzYMZSUlGDSpEnYuXMnZs6ciVOnTmHfvn3YtGkTjh8/joqKClx33XV49tlns3YtbGpqwuDBg3Hbbbfh1KlTuOGGG7Bx40YMGjRI+5cjhBCSX3gZsJ3OHV7dPeC8zV3TMa/0fETprc834oOFS3w8pT11i5rO33ft/bbnqXzPtkNVOSc6vBA4D0ocMA8KIYTkH36tBV7Y3NXv7yiEyLy2b58/Vrch61wrgTLpxe8POE9c9x/rnrW856QXv4+TB0tQ/H524OymJWu8fwETdgLFbtknbryM3xQoKSRXzX2EkPwmLIFy6/ONAJBZ0rHCzWfDz7OZxVBS8CJq3D7jFQoUkhrsnMYIIclDbq/z2r6Ngt2lA5Y5xLLFpeseAwBsu+VxpWubl2e8IvcdspUk6HWTiB//FgoURShQ0o1bo6JIISSZmNuusDDsOlRjKVaAfp+NE+PP2vqCBLG6OFkS5GezWtLJZ8LsQylQArC0rT/ZnHktUaURcGCMBtkMa7ceLOBvQkj8THrx+5g1tt3REiEEwdY9V2Hs9tNo2bkMAFA/ezW6agqzzjVmdHkSDZu7pmPXoZrM3+/M+aHj+VH4wsSJm1UoKQIlcB6UXEWugKJhNP9iKmbe1b9eaF6LjPMHJdnYlbXVb5TvXu6E5ALvzPmh66Av2m7z+1PRVVOI+tmrMyKltL0vS6TMGttueQ1Z5Mh8sHAJUKf2rPkuToDzZW3uM5PWV6bagmKuiN96bDFm3vWapRBJ2g9HshEmWTm3AdD/uwmrmGBe6V7Ma/u26yyKEOKMHDILOIfNCoRvyfDqbksryK3PN6L4/QuyRElpex8AoKumECfGn3X1SxH9dd2iJsdnEs/i5Dyb7wiREtUYxyUeDwiR8q3Hzu/1YzYfmh2KKFZyC7O1TBYvhJDzeBUc8vknxp/F8OruTDit+Kzdsrl5eVYswxTszt5M1pjRlRWiqyJQAOv2bZ6sAMmLsIkaMbHTne/FDi7xeEB4i4tUcifGn8Xcse1ZpsIPFvYfE+9967H+RqkyWyDxI3dUdWNjfBBCEkzboSrMvGs6mn8xVf0zpj5wadvt2Howe3nFTgD0WzDPWzHrxqJ/GWZO9nlL227HLtQA7/cLF3MuESvswoTnlVqcnHKiFCdeSb0FRSBmAuYlHnNFv3TdY5kGQoGSe5hniAB/R5IuRBsw13u5bxPoahterM/mpfd5bd/GyYMl+GDhEkx68fso2F1qm+BM1SoqnHaB85YdN4f7fMUu3DssuMQTgEvXPYbh1d30T8gzrISJDEUKyVeEn0VY4kMFVYHilKUV6A/5PXmwZMASj9fl2jQ4wrohltRmjW1H8y+mJlKgpH6Jx0xSTV1xkhbfGycRQwFDcpFJL34fwMAU66r12cnJ1MoSI/oKMSsXx4WzurzEKp8rJoTCqmHOwCqECXDOuVZTUrY0CxVRfv9Y9yywNuaHsYEWFIKlbbcPaOyiAfd3cO55A+LAzSqimySIFCsnxrZDVaFsakZyH2E9MdcNVTHutiQqCxhhfQaQ5ehq126sBEpcqKYjyEdBI5xkZcfmMOESD3Elij0XdBK1GLEiiQJFrMWLdfqoOhmS2wRpT071y+66cuSNLKYvXfdY5n2rnBxJTpCZL2Il6r6DSzzEklzKuhqFIBG5FQTmbJVJx5jRlfn/O3N+iEkvfh9tC5Nn6SLJQwxE5mWa+tmrAfhvC21r789qu+Lv4dXdmYFQpvj9C7IiJndVn08rny8CIOlsrtuAb+1enMiJDS0oKcK8AdauQzWxm1bt0CVQ5EYnOl9V5FTb4v+EJB25nnupt+JzSanrk178fmL7J0E+RP5s7pqeWCdZCpQUIA/2Ytad9IbvV6CYrSJ+ETPIJM4qZOxCRkn6cBLgZotI2PVF+LV967HsmfnSttuz0tCLjLFJEUV+yDVLj5icCmEVtf8al3iILRnntTnO54WN7iWcoMJE7iDrZ6/OXM+L1cWqkzV/Ppc7YpK7yO2jq6YQdYuaPAsD1bpcP3s1Pplfg12HajBSuq9IlDbvlvP7wDS/n50UTmR6lZPFmZefZNb+7Imsv+NwuFX1lYmLeW3fzmzUKMqn7hdNuPVgIwC1xHdxQQtKHhN2UjK3fS7cPhsUHdYS814fQXDqsM3n+TXDE2KH1yVMGbc6KLdXs7CRl4bMPixy/yAiAk8eLMmkxB+7/XTW/Ze23d4vbH4yDF01hZl72X23T+Z/OmCX5M1d0wek1Q+bukVNtsnjkoJYypGJw/JKCwqxRFdllDsrs0OcFWGsbetaytF1HcDfABGkTFTuR/GTDvyKEyEkRDsWgkClXasgW2qE9XZsex+AYXhn5xJg4flz62evxqGvXnVuRt/nKk6A7KVq4Q9itRt92LStvR+XrlPbIygu5pXuxa4Z/X6HVlmDkwgtKHlKmNYTO+uH7MFv9ua3Ord+9mrf0QI6hUWchC1QdN6PJBNVoapaX+Q26dTOnJY0zddQsRoGsQABwJTVb8a6E32Sl3msMPsHRQWdZImjiAjrulbhhVHh1JGqdrgqWC0JWW0J7/VediZz+W+d5IJQoQOwd6zqilXd0oFTHbKymrr5sAR9tkNfHYzh1d2YNbY9k8Jd7I58YvzZ0PN8yEktkyhW5rV9GwCyyoUCJQQoUNTxun266nWs8CNKwhQyZsHgdI7bear3M19DVaSEKUZU7itIilOvriUGv5j9KnJB0KkQlUBRuaeKwPGL3A7FBrDm9PthIXzzkiBS5OR3IlpHtC1jRlcs0ZwUKBFgNle6zVzyAZUN95KQ8TVpJH05yosFKN/qtEBuv1blYWUhjFO4hHVvv8LA7Vm8ihQdQurQVwdj7rW/yfwtb4wn0B3iXLeoCSfGnwWAzO7LceVJkTcDBJDlOCyeM47tMegkqxErM7O58UQ96yX6sbJ+6ELVguN0vqrIMS9nqSxvebXy+EH4GyVpuSZLTJ+LGJHLS/6/WXjXLWoCzoXrqqDje5uj5tz+drqOmaAiWrdgsrLsOfnR2B2TfVK27rkKu1CDAul4V02htme3Ktf+DM+IRaSI7y7+lS06mzK6JNn7d9GC4oC5woUZiho35jBCYOCz5qIFJUzhYUbMnLx6x6um3A+Smt/ps17Ej7AcWCHqi9VAKTtEOz1LVAJGZVIRZr1x+56iDOV/BWZLl5tzuvleOgWKCAUGgvmUeFnuUbXWdNUUYuZdr2XeFxviAef3BtKdRVVYJuZe+5uMxSIJSz1m4tzWhEs8ARGKOgyBYkUSRIvVd7UaOFR3QXU7N5cQDnZ2f5vPBbyJFK9+MF59a6ywEyjygGPGbpaqsiRkN1g5fYewxIrTd1BBpfy9XkfeERhAxqFThOm6OV6bfzeVcpVDi92Q+yg7gaPiz6RiAVEROn6WlOQyFMs/YW/3UbeoKeMDkwSSsN8aBYpP7KJQovIhUG10KoJGxUdGHFMVEiodcxKtKG44CY4wPmdGt6OuCmFHNql+TuW7hx0e74RqdFiQ64lyCLK0Z4fKM5qtYk4hwF6chr0IDjNukUde+0ArumoKQ4nukeuZvItzXCRBlMjQB0UDUYsTQH2TL7sOwq5BOjXU+tmrAWmAybUdfb1gJyi8iAwra4rXa5jRWeZOFhAZnfXayY/GToCUtvcpWyLiDDfu0tQ2nMpbHFP1GQp6b6t7WPUnDWMWYseRdVnvib5E/Dt0/9EB55hx8x/Ridd7FL9/ATYtWYO2Q2syUT46HWejTIaWNCGiAwqUhOI0S9Gd10DuwFRmtHbnJMkB0gq3ZRm387xeT+U6TgOf2bSvOkCqipQwcFo6sjtPlaBCJUhyQC9Lb0HQ/btZXc/8ntyHDN1/NOuYlUgxn9cwpj8d7KnaSwI9qxVOky4VEeHUV5rb1Oa6DUBdfzjupese0xbhIof6hkU+ihOASzwA7E2/ujsLr7MIJ9No2DMSL+ZmeQ0915Z3BLqWa4Jey2q262ZpCGt5wCtOOWCshEsQi4TZgdSPaPFSV1X8QPIFs0gBzosPq2Pmc1SdvIPi1cph119aZaC99flG30s/VvUqrH16clGYcInHA1GJE0BvOHJY4sRLojHz7ETnLDPqpSadTq06hI5duKvduVYiJeoB1SlM14yu31e0X1XriixmnCKSrMglYRIk8Z+VGHESLYIgTtt2WE3EgggTeTIl+qzmX0xFM6YOcGadeddraDtU5XvzQZ2TnrSSeoFiJikOsUDy86tYmUx1oStSQhdBHVnlZSSnTsvv8oyKGND1+yThN/Fi/TCfG8bmlUklSOivitXEijB92YJaTUrb+yzLwLzrsjnBG+rU7xlW1M7mrukA+vOa6LKcmIW9VTi7II72QoHikzBTkztdU7d/gdnHIcmzRN0dn5cZTtDlCBE+Ku4L2AuVsBPGBfmNw9gOIEyshIy4f9InADrQ4aAqW0rsrCqnai8JdUnc7+Do57sX7C4F7uoXJs2/mIp3PC7z/GPds9pyn8hCpG6slkvaJh+1ssDKv2kcmZNTLVD85Dkx/0BReaebCSomkjAL9kpSntnKuuM08NYtakIxzieISiN+LGJWyzVh7tKdr+jsn8zLOkKkDN1/1NFJNqhlMKzBUbRb89LfpBe/j5MHS/CBz7r1rccWZ/7v16ISln+JlbXEql3291XxTlpT6SRrtSU44E+gBGn8ZiuMF6uMToEiDx5hVkavM/ioBYkuZ1OnSBbVpR4/qEQPBfkuUeIkOryKlDBSu+tCtc35sX6ZLRxhYhYoTu3ISyi11ff2u8mgSoRZWPlQ/DjJhu0Aa7dMaszoGpCa/9bnGzF2+2ktItHL+J06Dx6rCiznZAh6LT+fl4WJ6jXD6GCjSAqlA/F7yS/xvt25Xq8t/9/r9zI/l/x+8fsXZF6CE+PPRmZZ8fJd4hAnbWvvdx0kxDmijFXOl891yv2RL5j9RobuP5p5hcGp2kuyrm2u/177FvHbyKJGvHRFCprr96Yla7C07XYt1wayv7PwH/FC2CnyvYqxONpL6iwoTruWuqE7/4h83biWiYDkzCbtcAuvNc+G3Cwhfi05usvJqwDQkRTOCqsEYTrFidflnaiXbHS0va6aQhgzugAAJw+WYHh1N0b+ZFjg6+pACAchIsK2pgTBnPjNLjpQ1FEvdcXtWiLr67y2b2tLf1+3qMnSIqH02QhCiK3GQ4F5TyFdMMzYArGnhUzcwiTMa6qSdHFiRh48nUJa7cJu/VjLwiojp7wyqugIZdSRRM2KpEViyQRpc5/M/xSzxrZj16EajPzJsAHfb3h1NwC9S6ZBriULkihEikpIsozTzFzeisNpR3kV7Lb1EL9f8fsX4FuPLe7f7XiO58vnJFYbwpa268ukG5S8t6DYhU2pNnidPifEHyrJymScQoJzNeGW0/f2KlLsfI50ROiY8XvNMKwouiwlwulRpEZ3IpfqmE7MIsWvQAkDK0vKzLtewz/WPRvKtgqXrnvM0348tz7fqC2LbRKhBeUcTmuVfsRJ2HhxHksTZnHilJrer5Usl/G6N5Cd5cnPEpKc9MrpPk7IM7awJgA6rysSexVYHDMLvimr38TWPVcBAMZuP63l/iqTK1UriVMiNF2hyeJ55GMqe/h4xbzZqx3y9xN1VETZhLHZ6QcLl2DSi93Kyzz5LE684kmgrF+/HuvXr8eHH34IALjiiivw/e9/Hw0NDQAAwzDw8MMP42c/+xm6urpw9dVX4yc/+QmuuOKKzDV6e3uxdOlS/PM//zNOnTqFG264AT/96U8xZswYbV/K734bMlGIE9E5JNEEnmSYnfE8ZutJlGVjlfTKa102DyRJMS2bEd/LmNGVsZrI1riZd72WOXfrnqtQ/P4FA4T0J/M/xcmDJYGFioro9rKEYxYpYYjEqP1eVCK9suva/bbn6eCdOT/EpetKYt/ZONfw1JuNGTMGjzzyCN544w288cYbuP766/EXf/EXOHDgAADg0UcfxZo1a/DEE0/g9ddfR3l5OWbOnIkTJ05krtHY2Iht27Zhy5YtePXVV3Hy5EncdNNNOHPmjN5vlhCS2uHmInbRMfJxL9fJF0QUUJBIICF0ZIEjomXMUTPimFy3RYSFKirROklj5l2vYXPdhowzrMCY0dW/5HPuNby6G8aMLhS/f0HGegIAs8a2Y3h1N6asfhOfzP80UROTlp3LQovwsbqXbusJgKwILSvqFjXFulfYBwuXDIjm8RPdkyY8WVDmzMn2HPr7v/97rF+/Hq+99houv/xyPP7443jooYcwd+5cAMBTTz2FsrIyPPPMM7jnnnvQ3d2NJ598Eps3b8aMGTMAAE8//TSqqqqwe/du3HjjjZq+VvCty3UlBpKvY/aYFn/nql9EHPh1dM1H7ASJ6lKNk6BxEg/ysaAWk6hwi5TzE0oumDW2HUC/74Bg7rW/QTOmYnh1d+Y4gIyZf17pXqCuP6FXWHVZjt6xwtwvDYX60lAuE6cjaPMvpmLr+KswvLobJw+WAADm3XI+iVsubv4XJr59UM6cOYN/+Zd/waeffopp06bh4MGD6OjowKxZszLnFBUV4dprr8XevXtxzz33oLW1FZ999lnWOZWVlaitrcXevXttBUpvby96e3szf/f09Lg+nwgbK37/At8dgB+R4ifZmmrnSB+VftJeDlY+OHJdPzF+sK+lHnENOxEhtwWnTTbl8M8wnA7jYtehGuw6VIN35vwQ9T9ZfS6a503sOlSDrXuuyvIdEN971th2xyyixowuoF1vOLJVDhQzp2ovSbXDf1ISENYtasLmrumYV7r3nHNsLI+VWDwLlH379mHatGn44x//iOHDh2Pbtm24/PLLsXdvfyMsKyvLOr+srAy/+93vAAAdHR0oLCxEaWnpgHM6Ojps77lq1So8/PDDys9oF06mA/NMzE7AyOfZJYez+r8dbvskpIUgOTrsyiyKstRxDyfnYFmIi/+7lZOViFHZDVjGyvonXyNJwsRPviHZQVu2hExZ/WZGeDT/YiqKTZ+z+97m5Fub6zbgWzXn06IHnVSpYmVdMSdbC4MorBYqdd9tKQiIpu4Kp9yt46+ib50FngVKTU0N3n77bRw/fhzPPfcc7rjjDuzZsydzvKAg27fdMIwB75lxO2f58uVYvPh8I+7p6UFVlXuWvbAaulNelLCSrvkJkU4TfpfIrHKqqIbLes1QG4T+zzsnrDPfyy4022ovDr/5VHLNxypI2gA5YZUIM35nzg/Rtlb9/sKEL4eeyo63w6u7YVQDXSgNnPfESWzYLePk8vKO/Fv6zcFTP3s1cO4z8v48YZIkEZ80PAuUwsJCfP7znwcATJkyBa+//jrWrl2LZcv6G35HRwcqKioy53d2dmasKuXl5ejr60NXV1eWFaWzsxPTp9s7CxUVFaGoqMjTc9o1bi/739h1vk6dWhRmU4qTflStULJ4UVlWc8oP4neGq2OPH6fPB+nkVD5rFW2WTx2rU53oshiw6sYexjsBdpftXw5akhEq8/Dt/gy06Leq3HqwEUE2arPbhycqR1g74tgRV2D+ja2eJUmOy0TDXjyGYaC3txfV1dUoLy9Hc3Nz5lhfXx/27NmTER+TJ0/GkCFDss45duwY9u/f7yhQdGHVMOz2v8m1mWGSCLuR+3WSla0lTp93Sl4mrAzy3iBW93E67nReENO0V3IxkkY3LTuXZV7AwDIJayM5gfBbEctHWZaUGV049FV/boJmUSJep2ovUbKwhImXPceC4tZe5GV4K7eANPvp6ED81n7L0VMm2QcffBANDQ2oqqrCiRMnsGXLFjzyyCPYuXMnZs6cidWrV2PVqlXYsGEDJkyYgJUrV6KlpQXt7e0oLu5fpf3ud7+Lf/u3f8PGjRsxcuRILF26FB9//DFaW1sxaNAgpefwkonOTrG7FZiqQGEFTjZu+/hEjRxh47Ss4mYRko+HKaad9urIV3Fjl0MjTN+ES9c9BiC7Xmxasgabu6Zj656rXHOn6I4CjCqaR3fdlXeqN5eFn2VccS4nrN6x89UMLZPs73//e8ybNw/Hjh1DSUkJJk2alBEnAPDAAw/g1KlT+N73vpdJ1LZr166MOAGApqYmDB48GLfddlsmUdvGjRuVxYlXrPxFxIzJj7gwLw3lK/kebujWgXnFq+DxK07Mx8OqhyoJ2KJao48TIUzCzp/xwcIlUtqB/vdEdMeu6hp8Mh+Omw/qXvbN1bYvB0h4XVaNa9KSBuT+5IX//V3lz+X9XjwCKzVnJ1D8+J7oIm4HWN3CxC0XQ9g4dTg6BYr5mm5758iYU/Wr1IGwBbJKhuN8FCd26dKjEmPypqYtO5dh0ovfB3Au2uex84EC5v2UVPuNpIfo66jXbmJSxV/N6vx8n5Tqxs51wsv4zbgmCTGYxrFs48UPISzE9xfr1Srnq5xnDl9UydOgC9XOKEiiLqfj5gyvcsZWIUzMgkUMOGYfFvF3nOIkiYOaDmQriZVvTlRiTM7QWz97NQp2l6Jgdylufb7RMrGeyqzfqi5FFcbslSj8U/z0rxQn8UCBcg6Rfll2mIuKJJgV3ToUKzHi5mwnCx7zNayiDOKOMDBj54yrmmq/tL1vQPp4EQHkZRYnBpfS9r7YrRZJqKthEnf5mhF1RfigmNPsy+fIqDhoB0Fla4Mg7VmHSPHjTG/3GfoaBsfPuJpqgeI3S6wu7Bq5jtmNqnXDCuHpbxeq6HR9L/d0sqroREVQ2DnNmWed5s942R/I6vNmAeP0HaISzip1MmkDeVCSFtFk9VsXv38BCnar50cx58OJg6AixaswqJ+92rI9qrRVu/PlaxP/+Pk9fae6z1Vk86lXa0muVNCgfiR2fiPmrJNWyNYQL5YRq8/pxC6ZmRVyxIaOzt3LNczPKKeNj5N8t5wkFbm/8psY0JwHSAduz6IzK61VJJlVm7Drn3UKtLQESfhFXh5uswlE+epcdYfz1AkUgZeEbVEh+x54aVQ6stdadShBHFzF9VQ/K+dosPuc1fOEYYES2STDmnn6ym4JdoppQ/69dWeo9lu3k7qpadT9NtvkQOycYoP8NqkTKHaFFYcwMUd8+G38QZ/dbjMxFb8Uq8/puLcVVu97mRWqlGvb2vs9z1a9oBpVZHaaNpuaxXtRdpBJsOSkCR2DoFx3gkYI+rW+6LSKRiGMkjBZzSXk3DPAwN8oSHnmvUDJhUrmlLnUD7qXSKw6GPGeldVFft+rT4rXZ7fL5CrjZc1eZ31R+S29dLhWe42IpUrO6PKfoLNRnYO7ynLpjiPr0DCmf3veXMqr5KWM2e7O41S//NbbvHWSjTKdsk7iMJ3K6a/tsIvgUV2KUX0OL8iK3ckBTkUohOEA7Ybdb+1XYOnCKsSW1pPocauTKpE04jxV7CYUVm3N6jk4UKcPp9886Dic0wLlq3ObLL983MIkaeGyqgihIosPO/Ei7+9hh2o5WN3DbY8bwD680gqnzjysHaidcMtoSQgwsPM3/+22xOk11FZladfpfuZ2tOPIOuV7eyHuPp6oEVSw5t0STxIqbpimTLt15LBNqFYRPUGEmJsTrdk/xy0k2A271PJexYmdJUkFv1EYZuLONkziRdTZ0vY+fDL/U8wa+ybeWPZnmeNuFji3+uNUp71YY0S7Gqr8CZIPWG0v45ecFyg6CyNuVJxk7Y7Fsb4bRKR4FVQqg7JT5ymLE/laYdUZswlcl6BguC+RmTW2HfNK96K5ZioA+3w+dpZGXfXSTfCEYVVOUn9vtyltWrD77uYJYFdNIYoP/FH5ujm9xCOTpMoq8NooOSvOxutyjVVyJvGvbOrWWc4qS1067+dmstfVDoTfCX1PkokYELbuuQrz2r6NE+PPYuZdr3nesdev4KVQpr+NCub+yGtfmNMWlO1b+zvPJIoTwJ9VI8hg5pT51e/13K4TVpp6tw5Q1fck30UfO8l0058CfxgKai7AruoaFJiOW4n2OMRFrkTwqOLmG0Sykcunf7PAFUqfy2mBAiRXnISJypKBrmyy5oicMH1d/OQziTJ8MYnOz2k3LZN+Stv7gPZhANTSuOsUKfk0CfDSluoWNYWa0DEfCNo35fQSj5eUuQIvDdNqF9AkoOqH4nc/Hj977KiEKuvA/FuIe1rtcxMlVt89ih2qKU7SS5BZvNVyp0pf4XVPm1wiH7c9SRr1s1cz1b0uzJEkYaHTYc1qycVryvmg9/ey944KuqJfvOIlTFpGNf1+GNYfCpZ0ITshmpP2CcwWV7v+zK0uCstLnBOBMBBtRuy/JVDxv7Jytmcb1EdqBIqO7cVldDs+6hZBqkLB76Z+cRFGh2iXKVcHfsM5vcJOMb3IqcStRIpTqnvVCDPdWWidrqeSGl11408V5A3uvJAP4ixsUh9mrEKQyhzGZnRW+O04BG7ZW512KFYdjFVm/F5T22djbYnwamFS9dFxOmYuF6+hkuy8SFyEvfygI2ze7hqqGyTqnNTJ11GxmnB5JzryTqA4iRGnSm2eZcTtd+K1A1AZQN024YvCCVR+BierxanaSwb8DjqWwrxsjKgrf0OQXaGdMOcAohWFOGH2OzFvJCgf04WX61rVXydrig6RZKZuUZPrJpzmY+ZnY1s8T9A8ZTntJGuF3R4R4m877BwazeutYc2MW3YuC1yp41x+sXLIlVPlq4gg+VxxPT9l7rasYvVcbnhJ6e90/zBFIGd26cZL/+HmvB00OMDp8+Y2Lfd9UdVhnRNQihE1/JRT3llQgGxLyYnxZzG8uhsFu0t9XUf+Nwx0ZsIduv9o1t4XYidRp2Ufcaxl57LM+eIzdtE8Xv01xHeUr6+KuF+YDsuyhcOL46oXQRiFMyw7SuIX1T4ujL7QygIhY1evnawzureDcNu1mG0vHPJSoAQhjAZoJULczIRB7iFQXaIQn9Wdktpu2UQ+riJ25PN0/z7mZ3BbBrP6jHhPPjdMSwnFCbHDi9OrisXYL7rFgeqeWVaZclWTOgIDo3Ks+mm2Nf+07FyGa2Y8rHx+3i3xCERFK37/Al/WEyv8zN6jqMxuyyZOx+xMq34HWKvdSxvGLLS0nnixxKiklPeDm1OxH+xElo5nZ+dIoiYJuaBkwSDagNvyr9ccRCqRRSRa8tKC4qUxta29Pyv+XXVTOh2hbeZK7yfayK7hNIxZCITgk+K2VOG0jPPR16sz/9fh6Opm1fB6TbPlwy6XTNyh1oT4wSmHidlXzyrnSdQ4+aXYHbNy+g0irpjfJF7ywoLiN1unCCkzb4qmuhmdG2Gb4q1m5CIXgsogag7p05V11ureKs/jdeA3WyWC7Kxs9XkVR9oo09871R06yBJzojEZFWdYK8ESt0gB/LUxp81Brf5uW3u/p8gdEg05L1CCpBJ3atDimn7Vt1tlt0pTrZLtMSy8OLA6iRL5X/N7dimy5XI2+4Oo+Ie4iRQvETh2qPimBLmOGyodpBCnJJ3I1gM/iSVloWIe3OMSKQ1jFlpGADrVdS/Rm/L3snLO1RFdSfyT8wJFNQxVWEnM57qJlChRbQhO/gx+GpPX6Bor0aAjbbvoWD76erWrUHGybMj+KnbWERXM17G7RlCRExSrKIgk1WsSDV53AHc6HvYeUm6IqEJzu3cT4C07lyklWwOQyXci+kyzRZnET077oBz/fCEGFRVmiZSWncuyOmerymoXsirONXfuXvefCGL2b1t7f1bqZSu/FF0DnWjsQz18xikjrRXC78Sq85PXuVW8763wkmROxz5DceC1s4zboZHEg46omzC39NCFk0+KOC6HRVh9BycRQ3HiD7ffxQ8FhmEYga4QAz09PSgpKcEV96zEoKILAZxXwyqIgnRKcWw1Aw3quOoFq2e0eg67e6ma+oMM2iqb5HlFFiul7X3K2XGj9AVRjT4yP5uqX5AXxzyr37mrplB5FknyBx3Le1bOpap5RnRufhl24jaKEL14+Z1On/4j/s+/r0B3dzdGjBjheG5OL/H8f7/ts1X+ZuoWNTkm27HaydKpkw97lipMj3adg8CuYqg0QB3ixK0cvAyUskVFVQyGEXocFKvlJy9Oy17Xvs2maYoTooqXjK/mYzLyUqqu9ug06NE/JFmI38HsDxg0PUROL/HIeHECs8Mq9FcOQ9bh0S72ehCoNDDz8oj5GbwmD9I101FZpzYLP7eZmXgvaaJDxmveGa8RVV6wWj8nJErMSQp1oDojp1N48lCN7FQhLwSKyvKO+Rzxf/MA6jTYB8myaPYrkTMWelmm8XJfp4ynKgLALSOqXdSR2/u5IEIIyTXC2vDPDS+Dj45JnmpWWafPk/gY+u4x5XNzXqB4qWx2FhJxTGDe0dJPgxrgw+KQhE2c65YcyakDshM6ujOiysfkjsmcOtvs+Co/c64IEy8p68NOb2+F2/4ghKigK729G7qWxf2IE7aR+DDvcXbq8gpAUaPktEDZvlXvWrudALDzbHeyfAgvcq+N0qmDsDoW12Av+6B4XafOFZz217EKrzZvvKg6s/TbeYq6xf1BCOBPXERpbZH9BaO8L9tGdNhN1LxuxCrIaYGiE/OSj1OltoqScFP1fn1k3D4XdRSLlXOsnYBTXUZSPTcu3BqVnCbfai8iJ4J0nsI/imHFBAi+9OEVv0JDhzhJQoZbcp6GMQsd01X49XvM6SieqDAneLNqGFYDjcqasNvgojr4eB0Y/XCq9hJLr/8T48/6vqZOr385gZtbinovOGWnFfeRxYnYamDo/qOO0QZBZ3b1s1cP2H2VECAcwZ+ETQNlkvQsxJ4gdTE1AkU1JM1OTFhtu62Kaqpl3egULUKcWDH32t/AmNGV+Vt0ZFEmNvObit5tt2eVc1t2LrMsa/N7YQ0a4hkIiZKoLBhOYdBusF0kA79jQWoEig7MeQFkkeImWNxESlRbfQexMNg9565DNdhct0HXI3pG/h5W5RzEmmL3WafrWYlhK5+UoFYPziCJGVHv3NqECmZBLbf/KOueXToCt2egOMl9cjqTrEomOjfkjK1uzq/y+QKz34FdBljV2YbftVW7Z5b32bEaVIPsIyN/R2NGFwp2lw44p/K5g0rXd7qfSjZZc+SQl0RvVtdzWjM1Z4YVf7tZrHSmgpbrGROzETOZbSzO1eOgmV/jCmF2wu2ZKFCixWoSZlXHvGSSpZPsOVTEiZkgS0a6zgf6G2rdoiZfA5VqeKz5HHOlsxIn8nlBhJCKsPK766rdZodedk6WzxWC0EqsWKXv9hOBI+fUIcQr5joc1x5TQff9icrqTPwTtG55WuJZtWoVrrzyShQXF2P06NG45ZZb0N7ennXOnXfeiYKCgqzX1KlTs87p7e3FggULMGrUKAwbNgw333wzjhw5EuiLWFG3qCkTlSO257bbUlsX5twpXgcSP6mBhenV6vvtOLIutA7IqoM4Mf5sltOs+P5OSyVWS05uDnlmC4tdman+tkHKSNzb667QXhH1meKEeGHo/qOOSyLmtmOVotzL9hOqyJOKJFlmiD92HFmHlp3LbMcu0ad7SQ/iaYln9uzZ+MY3voErr7wSp0+fxkMPPYR9+/bh3XffxbBhwwD0C5Tf//732LDhvE9CYWEhRo4cmfn7u9/9Ll588UVs3LgRF198MZYsWYJPPvkEra2tGDRokOtzyEs8N9+2fsBxvxvs6cS84ZuVMDKf52alUMFpicguJbUXMSQsLmKXYjPGjC7MGtuOXYdqcPJgCYrf79fAXjugQ1/tN+6N3X46856dSdfLZnxBUSkrr0s9KvXRavNKAZd4iBX1s1dntQ3RN7htbgn474tUl4LcElLqECy0osSH04anXlw0PC3x7Ny5M+vvDRs2YPTo0WhtbcVXvvKVzPtFRUUoLy+3vEZ3dzeefPJJbN68GTNmzAAAPP3006iqqsLu3btx4403Kj/PV+c2YfDgC718BXTVFEaS2Mp8fbv76Y7sUBEn4m+VTserg+jJgyWYV7cX80r3YvPY6diKqzIiBbBPiW9G/oxAWGZKs412sZmo40LFQZCQgcuR6vl8ZLy0L1VhoZI1O4hIoTiJFzEpD/o7BIri6e7uBoAs6wgAtLS0YPTo0bjsssvwne98B52dnZljra2t+OyzzzBr1qzMe5WVlaitrcXevXuDPE7eEVS8qC4X+Q3RtaL4/Qswr+3bmNf2bew6VINttzwOY0bXAIdat9wpVmbf4vcvsBQuKugc0J3KRWdot1jSMVtPVEUeSTdWdVFlo0uv7d5vfhTznl0686wwL1C8OJX/17/wgPJ1fDvJGoaBxYsX45prrkFtbW3m/YaGBvzlX/4lxo0bh4MHD+J//I//geuvvx6tra0oKipCR0cHCgsLUVqa7VRZVlaGjo4Oy3v19vait7c383dPT4/jszmlqk9ah+7mpOp3+cKLuNFtxZEdZjePnd6/1HPu79L2PhwaX+Lrul5/OzmyR+fvbldeKuLEbXnHaSnHDJd2iBe87MKt6jxvzijtt53p7pdpQUkGVkLldBR78dx3331455138Oqrr2a9f/vtt2f+X1tbiylTpmDcuHH41a9+hblz59pezzAMFBQUWB5btWoVHn74Yb+POqDy50Ll9ZMC3qvQCCJMrDoUq9nP1j39Szz9FpP+Dmzs9tMDZktOmwsGIcqU2A1jFrqKFKe6ZxYnIpW9+T1CvOLmc+Y1mkdX9E8Y7TMX+neihi97+YIFC/DCCy/g17/+NcaMGeN4bkVFBcaNG4f33nsPAFBeXo6+vj50dXVlndfZ2YmysjLLayxfvhzd3d2Z1+HDh7OOW4V+6kyhHhZeoj68OrImBbEkI/61c2AWHZXV0k4Qs2+UiaVO1V7i27RsZTkxh45TnBCv7DiyLhNdIVDtH9wSHIrooCBQnOQfOqMZPQkUwzBw3333YevWrXj55ZdRXW0dySHz8ccf4/Dhw6ioqAAATJ48GUOGDEFzc3PmnGPHjmH//v2YPn265TWKioowYsSIrJfASpBYhczJ/w87HFQXSRdYZtw6G1l8mMWIk5AQxw59dXAmukdGJRW27o5QpwiUxYkQIWZhQnFCgqI6eKv6rYk2oDNM2OyLojqx0J0ugninYczCzD5kuvAkUObPn4+nn34azzzzDIqLi9HR0YGOjg6cOnUKAHDy5EksXboU//f//l98+OGHaGlpwZw5czBq1CjceuutAICSkhLcfffdWLJkCf793/8db731Fv7mb/4GEydOzET1qDL0XcWFLAwc7JMiUuw85u3MsCrryDqEjViD9pOTJQhukT1+nWRlglhT5I7ZnL9Ft6CkKCG6cRrErSL93M4R6LZQ+k28SPILTz4o69f35xypr6/Pen/Dhg248847MWjQIOzbtw+bNm3C8ePHUVFRgeuuuw7PPvssiouLM+c3NTVh8ODBuO2223Dq1CnccMMN2Lhxo1IOFJ2YRYqT/4CVoAkSsWEnkIJmddU1WKp+3vq8/mcQndaJ8WdthYWXtWyz17/V+14Jms1SIH+HoOF1FCYkbOycYGWn/DgsuDq3+SDRIcYz3ZmJPQkUt5xuQ4cOxUsvveR6nQsvvBA//vGP8eMf/9jL7UPHzsmxfvZqwKLBqjhFRkEU6aqtoons9rExYxYnsmOcLkc7QK1zMzvlmZ1z49hzxMoZlpCkYO5f7BI+xgHFSbJQ2SLk9Fn1vpW7GQekYczCzMvpHKe/vWLlvGblKGz3uaBp3eXrOT2TnbNrXGZbeV3bzSTtlJdBd84GQqLEKR25wOzPZ7fUG2YqfJIbuI1nQSxx3CzQhFzYwjrSsnOZkqjwIlL8YrfTrlksqDi5+a04bjv9ysesLBTmpRo/+U3M6MpgqXo/u2N+s7xyWYfExY4j63z1T247nIctWGg9STY6lghpQXFANNqkONTKmC0Zbg6tOkyxKmnvVe4TpBMzD/4z73otk6lWRRhYCSSddNUUUmyQnEPXUrXs3BqWlZGWy+TiJ7DCCQoUF5IiTtyie+QOJiznNpWMt1ZY5Tfxukzidp55F2Wnz4fdwdGfhCSdsH3nvC7nyn2BCBm2Cx3mUlJ64BJPHiA6GytTrSxsdDq0uUUA2OF3CcQ8MxM0/2IqCgDgnDiJMnMsIfmA7kmYWxuXRYfKzt5WO7+TdEALSo4Qd9K2oI618vNbWVN0orrUQ1MxSTNmK0oUEYlWu7y37FymlIzQ/FkKlvyHAiVHURUMbt73Ou9tFyFktzwlRMqJ8WfxyfxPlUWDlYVEVxI3QtKEX1Hi1v9YRfXocGqlY2xycIte1QF79DzAKvLIiiAOTEEtOLJIMVtTzFj5p8jvqUTx6ErARkhaCDOvE4VFfhGVbyYFSgyIDbyC4CQYwjbV6hIrgrHbT2PkT4ZZnqtiVfESYhwFjOIh5DwUJ/lFlIEjFCgxolNIeEnbrxMdyd8EOlPWh3EPQvIdFSuKl+VlQoJAgRIDYSlQs9OYlRNcFMLFatMxqx2mrTowL9YQ8fLq7BqmSKH1hOQaVn2Crr4iDJEiX5OOstESddoNhhnnKKqZYKNc7rHar8NttmUX/uwlf4L5M2bBYvUsfsOdCclH5BQFunxRwtqnZ8eRdaifvZrtN2LiyAlGC0rM6BQQcZtU/UYLyRaVofuPahEn8jXcNjgkhES3LKwTWizzGwqUHMBqJmI30DqFfon3o1rq8eub4lVEmHcnthI4Uey6ys6SpIEkiPyWncvoT5YCuMQTE7IZ1e9mXV7vZ3Vfr593Q3fn5ZYyOymdlCxO6hY1UayQvED0EXL79yL2w4zgYXRQdMS15QstKDlA0DT1OiqXH4uLX7Fi/q5exImXdekwtojnPjwkl8nFZR4SLnHuR0cLSkJws6J4HeydruXXCc7qM26VNw5zsNN+PHYOuUEc7mgtIfmE3BdZWXlVHeAJCQotKDESV+4ScW8/ylh8Lsw0x15ETZBQZR1QnJB8xC70mJAooUBJOV5FRlSd1ND9Rx2tGuKYPItT3cnYnDbf7w7ISfF/ISQO4nKWrZ+9OvMi4RLn8g5AgRI7bhVAV5ZWnViJFN3PKK4XZMOxoMng7BA7r9o9C60qJB9QsaKI9hWlWJHbHYVKfkMflIRiHvDtdgROAmEJKFlIqAgTL8LDKqmbmy+KivCgOCFpI6gTvx9EfyDEifiXkT36iNt6AlCgJAKz02qUjT2J6azF93frbOyWZuyy25o/C1iLIDkSh4KDkGxkh1mR0dpOpNTPXu1JNMjXVfmclVChSMkfuMSTQyTJemK1z48uhu4/GmsnI8KPKU4IcUfHhEp2uvdzPdFfdNUUom5RE8P9NZAEp2gKlISx48g6WyGSRF+UMLLSun1P0RkJ64eO8pIFUcvOZZyFEaKImx9KVD4i5jZLkZL7UKAkECeRovMeSVDIVqh8d7N1w/wZLz47VmKEjncEYD1wQ86JEnQCpXsCxii73Ic+KCRR6NxFVayL2zm+WnVgYpfU0sBPQfIF+jWEi980+lY4LcuaLSpcwk0+tKAklKRaN8LA7+xLDBp2n1e5ptXAw5kXEXC5zxqr/klYU6wmBE4RIVbtN4oyZ4hy8qFASQhWDdhqeUKHGTQp4sfcMel4LqcwYbPw6KopDKWD4to3SRMikkdgJ/DN2afjFAh21lOSLAoMwzDifgiv9PT0oKSkBDMq7sHgC/zvoZJEzIO0VaNx8qswdxZu19eJaty8VfbXIDMmUUZuOx+b0XFvK4RAoQmZpIWGMQtdJ09OEy5zmHJYFhSnnca5lGeN7nwop8/2Yfexf0J3dzdGjBjheC4tKAkhioyxSbGciDT2QTboM1/PCbf76GyAnIWRNCL6FtV2bdXfiXYcpkhwmjRQnAwk7mRtFCgJQTROc4UwD75u1pNcQH5Op5TxKsiCQOQvEYiysssQK849VXuJNmGhS3QRkmu07Fw2oA3KOPVPSdzSg8QPBUrCMDdSs9XDrhHbvR9WrhKvmDsgp44sCLLYUe3wdFpzCEkzwqlYnkiZ2xaFCFGFAiVHcBpEncRJlOw4si4jRFRmRLpMqlap8c1lpSKGgjq30jmWkH52HFlnmVAx7okSyS2YByWBmPfmccLOKTaujsBuTw75GXV2VHYiR+4UvczYdDnK0UGWEOv22bJzmeOSql1gAMVN+qAFJYG4DahmC0UY4bpBcNusL8rncytLs6UlyFIPRQkhash756gSt8Nm2ghS3rqW8ShQEkr97NWeKoiYYcQtTsQyT9yIHVbdCMsXhhDijmrbS0KfQtTRtVULBUqCsNtPxuxbYTXriFuYOOHFL0UHLTuX+b5PULFCKwohanhdSqVISR7ybxJG/04fFPQP7kk1H5r9UeQBlHH7aoiEbE7otKJQpBCihtyHMYdQslAZF736+HmFFhR4c0qNg4YxC7OWIrg/iBrCgTjqJRxG8xDiHXO/xtD/3EPX0o6AAuUcSbCgiJ13rRqmaLwUJt6xKjOrhsR8KITEjxyeLNqj7oGP6CPM3ybvBQrXLaPHSURF1dGYnWTlqIGumsIB9UKnMLFa4qH5mhB1zDlUdPTj5s0KiT/MPoXyfkpW5wXBk0BZtWoVrrzyShQXF2P06NG45ZZb0N7ennWOYRhYsWIFKisrMXToUNTX1+PAgQNZ5/T29mLBggUYNWoUhg0bhptvvhlHjhzx/PCnLq+wP+axcOIUMk7P6vW5kjQQJsEaIS/duaXi1o281JOEsiAkF9ExqamfvZqTVR/I/afTOCV8UeTUFzrw5CS7Z88ezJ8/H1deeSVOnz6Nhx56CLNmzcK7776LYcOGAQAeffRRrFmzBhs3bsRll12GH/3oR5g5cyba29tRXFwMAGhsbMSLL76ILVu24OKLL8aSJUtw0003obW1FYMGDfL0BawSlakUTlIqq+qzWvnJmMVIV00hSrU+nX9adi5LpC+GSBJldu6S/Xt0k8RyICTpiLYYxOqRaXvnJggtO9cN6De5bO6M1yAStz2XTp/+I3BM7VqeBMrOnTuz/t6wYQNGjx6N1tZWfOUrX4FhGHj88cfx0EMPYe7cuQCAp556CmVlZXjmmWdwzz33oLu7G08++SQ2b96MGTNmAACefvppVFVVYffu3bjxxhu9PNIArApHiBj5X5k41zdVvaCFSBH/tyJp+TzMzxNlOX/09WpUPnfQ8ljYnucCeamHkT2E+MNvW7WbGFhF9Ylz5XYqhAwFTPZvYFV+KslF/RAozLi7uxsAMHLkSADAwYMH0dHRgVmzZmXOKSoqwrXXXou9e/finnvuQWtrKz777LOscyorK1FbW4u9e/daCpTe3l709vZm/u7p6ck67tViYh6gZCuMXer4sDA/u5PISIrVRxWxrCF/p6iipUrb+2zvJd6XZ1JhdUJ1i5ooTgjRgJdlUlWrZVdNYda5dYuazjvnclk2g9v2BE4EGbd8O8kahoHFixfjmmuuQW1tLQCgo6MDAFBWVpZ1bllZWeZYR0cHCgsLUVpaanuOmVWrVqGkpCTzqqqq8vvYSoQtBJySlpkjSRhZQghJMy07l2X6QB0+dqqWZnFe3aImLtMi+3ewGpPCGKt8C5T77rsP77zzDv75n/95wLGCgoKsvw3DGPCeGadzli9fju7u7szr8OHDfh87Cz++KzqQf0jxf/k9oeBzXZjIDqlxbPal0pmFHbotW0/Y0RHiD9GOrCweXtpU3aKmAf2tjOivzP+SePC1xLNgwQK88MILeOWVVzBmzJjM++Xl5QD6rSQVFecjbDo7OzNWlfLycvT19aGrqyvLitLZ2Ynp06db3q+oqAhFRUV+HtUSK0/jqMRJ2PdJ0nqpbBYU4qRhzELLFP5h4bQ7cZRlRWFCiD7EUkyp9HdYiAkjl2v7haLo080CLwwx58mCYhgG7rvvPmzduhUvv/wyqqurs45XV1ejvLwczc3Nmff6+vqwZ8+ejPiYPHkyhgwZknXOsWPHsH//fluBopO4w4mB8z+knYXEznxGcgsxu6M4IUQfcv/ZVVOIlp3LtAkHt36WbbmfqMYjTxaU+fPn45lnnsG//uu/ori4OOMzUlJSgqFDh6KgoACNjY1YuXIlJkyYgAkTJmDlypW46KKL8M1vfjNz7t13340lS5bg4osvxsiRI7F06VJMnDgxE9WTBHQ7y3oVRrIaVa0MSbKeCMzPlGuOvoSQZGBeLo0SOastLSnWqOx55hVPFpT169eju7sb9fX1qKioyLyeffbZzDkPPPAAGhsb8b3vfQ9TpkzB0aNHsWvXrkwOFABoamrCLbfcgttuuw1f/vKXcdFFF+HFF1/0nAMlCDoVYBi7ONJiEowkiLW2tfdnXoQQfcTdptJsSREJL63QPW4VGIZhaL1iBPT09KCkpARfvmEFBg++0PPnrcJfrVCxoJjDl1XOk5/BDVVFmoQBWYUoQnut7hd3+XDWRUj46BIOssXE6m8gfpEUN6Ks5bJRGa9On/4j/s+/r0B3dzdGjBjheG7e7cXj5KFtjpQxvyeQc6KY8WopsdunwI0oU7LnO2bPf134iSAghCQflQkk2/N5nMYrc6Tq8c+rW1lyXqDIu/yanaXMobxWWBWsWYSYo37s8phYRQc5iR078lmYCEtGVBYN+T46OxRzqCMhJH+xGz/S3Pbl5euw0mIEyiQbN9u3upvYghSabPVQFRhumyWJLLZmEeLFRJbrhLXfjYrJNa6llrSbgwnJV6xESlrbu9UKRRByWqDY0bb2fsv1MS+EtV+L3X5AQYjbvyIpyB1FWB1EmmdMhBA1rPb2yVfE+GPuG63GpfrZq/HC/56Pkn96UOnaOb/E40ZXTSHa1t6fWcpR9e3wKyBUxJCTX4pXMUVx4h2/IoPihBBCrBFirG3t/dqSY+atQDGHeOoeyM3p6uV/5XPcECLFLu/J0P1HBwgZq/fSjnmm4iQmdM9q0jBLIoQQN3T3hXkrUKwQmx3pdOaxEide72EWKU5Ou7IwofUkGzuRIixoAj8bjtkJHooTQggJh1QJFJ3ocGa1WkZKg5NslAhhYd6I0c81vB4jhERHEttiEp8pl0i1QAlqSZHFhFsMuB9ka4l5WSfqnYFzBTuLhvA90ilOgsLOixBC7EmVQDGb9nVaK/yKEKc0+XbOtBQnzphFipxNNikWKooTQvSR5PaU5GdLOqkSKGHkGlG5lpf7CTFi5whLceIdKyuWkx+K112I2QERQpxgH+GPVAkUXclj/NxTxm3Zxy5Ch+IkGHWLmkLZcdN8DxVxQ+daQvTAwT9/SZVASQp+9tmhOPGGkwAQ2yJY4SfCx+5e7DgJIcQ/Ob2bscpuiDLy4KN7ucdprx+BGBTrZ6+2vT+tJ3rxEx4s/z4q0BpCSDzk2iSAfYW38Tu1FpQodgu2u77X/CUUJ7mB192Nna6Rax0vIX6pn73al+UyF2G79kZqBEocDUD2NZFFid2zMDusfuKYsfgVGebz2ZmRtOAn6WSutg9OQtRJjUARDSCMbaH9Xs/OwiKHHtN6Eg5OnYNdOLLVe04CiOZcQtxJc0ZsihRnUiNQzIQhVFQQ1hMx2MlhxWbC2E2ZnMfNtOyWiE+3OHH7jDzzYsdGSH7AtmxPqgRKXEq9q6bQdjAUIsScsM0p0oR4w27gVxGocmi6F0Hrt9ORN7h0ux6FCiH5AduyNakSKMDAmXAcVhQSPU4O0VbC0UocOl1DZ+ciX6tuUZOr/xQ7N5JG8nEJVUxk5VeaSZ1AsRp4dO7Jo3pcZIo1L+1wWSd67CxV4v3S9r4BnaHdbx5Gp6laPylSCMk/0ixSUidQAPvsrlFgl8Jevj+XdvRj5fha2t7nWtbyUpt8DfPnVJdlvOBH7NCaQtJCmup5WkVKKgWK3WAVBJFXxZxfRf4/w4jjRcdmgS07l4VuWpb9UPzeK02dN0knSdn40y+5/vxRkEqBAui1UlhZRbwmgmNljYYodjQWAkOXkPF7LTeRQhFDSDyopjGQSaNPSmoFCpBdIeIKOzbD5Z3wiUKkhIEuwcNlIJIPxGGR1tFvmC3sVpNZ83eT/06TUBkc9wPETRw7HJP4CUMIRhFV0Lb2fk/iQpwrns38NyG5TFw71Id1PzvRldbgidQLlCDIlelU7SUYuv+o54pkPr9+9mpaUSJGDNpmB9h8/C0oTEi+sOPIusgtCZzIRgsFigbSqm7zja6aQpS292V1evL/802sEJLrDN1/FB99vTruxwiMleWEKShS7oMiDzhuTq0qOUvkv4VPi3wPp/NJfAirQi76pcicGH8WJ8afjfsxCImMHUfWZdptVO1X930Y3WkPLSgKyBXIahlHfk/8X1RiOxMkxUmyMC99JNUJrX72apTC2dRs9jchJA3kgh+KmzOsGbtxIi1LTakXKHZ+I04Vx6/ipVLOHayWc3T4pDSMWRhoh2r5/mYBUvy+vUG0blET/U9IXtKyc1lmQpEkkSI78Oa6dTYuUr3EY0dQITF0/1HLgUxsCEjrSW5i/k29hPuJc3X+9laCQ37PfJwWFUL0YhcmLAhLnKRF8KTegqILs6hpGLMQQ8/9n4IkP7ESLHaWF7vPBMXNKuI1LJmQXMWr9cRsAfGzbBOXdSQtDvu0oJxDCIwwlmFUr8kloNzHyaKShE6FYoXkI37alhAX4uUmTsyWEnG++XPyeWmxdIRF6gWK8AcQVo44rR3i3g1jFsb2DMQ/Tp1knOKEvickDTiJAqclGD/o8nVxGm/sjiVhohMVXOKBvaNsXM+RhGch/rDqPJLQoVCkkHxHdpY14yQogogNc9r6sJ10k9CXREnqLSgyVrlOdF/b7vpc3iGEED34XVpxssKoiA/V+7q5FJgnqUnZKy5qaEE5R1wCQSWDIEkXYonPSzhyw5iFOFV7SepmWIToRPYrkS0iZuHhVYiYEVujEGdoQYkJO2sKKy0R0BeJEP+E4aCqq392uo6V9QRI5zItLSgxoVNBi3VXzp5zFzsx4jVEOR83OCREFbMfShC/EKsQYi/+gZxsBsezBeWVV17BnDlzUFlZiYKCAjz//PNZx++8804UFBRkvaZOnZp1Tm9vLxYsWIBRo0Zh2LBhuPnmm3HkyJFAX8QvccxSw6q4nHHnJk6/Gzs5QrxhFuhuydTM51r9P05K2/tSaT0BfAiUTz/9FHV1dXjiiSdsz5k9ezaOHTuWeW3fvj3reGNjI7Zt24YtW7bg1VdfxcmTJ3HTTTfhzJkz3r9BAOIc0N0GHtVnk2cLjP7JT+wyE1udRwixx0106HJE9bKEY/VZ8fm0OscKPC/xNDQ0oKGhwfGcoqIilJeXWx7r7u7Gk08+ic2bN2PGjBkAgKeffhpVVVXYvXs3brzxRq+PREhO4ccJVgU63hHijuqyj52DrBmrXe51tcWkWHHiIhQn2ZaWFowePRqXXXYZvvOd76CzszNzrLW1FZ999hlmzZqVea+yshK1tbXYu3ev5fV6e3vR09OT9SIkFwnTaic6xIYxC7ncR1KNmzhQyfSquiwkkIVJ0M1maQ3vR7tAaWhowC9/+Uu8/PLLeOyxx/D666/j+uuvR29vLwCgo6MDhYWFKC0tzfpcWVkZOjo6LK+5atUqlJSUZF5VVVW6HztWgmSx5Yw5dzCLBgoJQsJhx5F1Sn2j3V46Xi0Xou9mf6wX7QLl9ttvx9e+9jXU1tZizpw52LFjB/7rv/4Lv/rVrxw/ZxgGCgoKLI8tX74c3d3dmdfhw4d1P3bsmMWJm/mfg1u64G9NiDfc+lCnPXPi9vug0Okn9DwoFRUVGDduHN577z0AQHl5Ofr6+tDV1ZV1XmdnJ8rKyiyvUVRUhBEjRmS9gpKkDt8ty6yMkzBheGly0VHf7K6RpLpMSJJo2bks1AzhqnDJxh+hC5SPP/4Yhw8fRkVFBQBg8uTJGDJkCJqbmzPnHDt2DPv378f06dPDfpycxmogYsVPPjoFBMUIId5QdUb3akWRhY8OvxOBeT+2NE88PUfxnDx5Er/97W8zfx88eBBvv/02Ro4ciZEjR2LFihX4+te/joqKCnz44Yd48MEHMWrUKNx6660AgJKSEtx9991YsmQJLr74YowcORJLly7FxIkTM1E9aUZ3ZAeJj7DEhEoU0KnaS5i0jZBziLYitoSIA5XIHk44s/EsUN544w1cd911mb8XL14MALjjjjuwfv167Nu3D5s2bcLx48dRUVGB6667Ds8++yyKi4szn2lqasLgwYNx22234dSpU7jhhhuwceNGDBo0SMNXyl0oTogX8s2aUreoCUA6U3qTZCCsKHbOs05YOcp6CTemOBmIZ4FSX18PwzBsj7/00kuu17jwwgvx4x//GD/+8Y+93j7vEBXYjzhhhSZ2DN1/lPWDEBM7jqxTsqJYiRMhNOxEh9t7cfvB5CKp3SwwCdYK0UhUn8W8NilDUz7JddrW3k/rCQmdoH0/hUZ0pFagJAHVFOby+YCzUCHJIAnLL0P3H03EcxCS6+gUJV7677RPPFMtUOK2osR9fxIOSRMFzJlDSDay4NCV8yTIpJETTmtSLVByDQqaZEMhQEhu0VVTiBPjz2oRKap+X2ZLOMWJPRQoESJXSD+VkoNfMskVYZILz0hIlKhG6jj11162KlHt++POZJsUPEfx5BvCqzsshMd3mCo57euUcZBrgz1naYRkI/rlsdtP+76G333U3MaE0vY+9uugBSV0rCqiH4crMSDSg5x4RQ6LrJ+9OuanISS/UbWkEHdSb0GJArf4eZIb5JrVRMA6R4heRF9uZwnh7sZ6oAUF0TmfMmkPiRPRoeaq0CIkTvws5/j1OeTyTj+0oMRIw5iFnsSRVSVnRSaEEH949Q+Uz/WzTON2P/bn2VCg5AiqKZoJcUKYpr2KY0JIMJgF3Dtc4omRoAMEQ9GiI1+WRbi0SNKOXVuOo21QnDhDC0oI6IjasUM2EXrdbZMQGdFR05JC0krUzqxil2QKEzVoQTlH0E7ayRlKdq7yukGg1T1I9HAQJyS/iaJv5aTSG7SgaELF2Wro/qNaBjpd1yHpwyrUnf4oJO1YWb3DEiy0nqhDgRIhQSs8KzYJgmzJk3PzEEKCI3wCnawk7MO9wSUejTitYzJJG4kTsxChMCHEHraPZECBohFWapJLsL4SQpIMBUrI0GpCkoDIIksI6UdVoKu2m9L2PjrBaoYChRBF6EhKSLqgqI8XChSJoAOQVWWWVToHOJJUuMsxSROq1hMug8YLBYoJigjiRC7WDz+bnBGSFtysJGw38UGBEhEMLyNJhvWTpJ2w/bTYxrxDgUJISuB6OiHWFhHmBUomFCgW7DiyTmtaeSpnkiQoVEhaEftPiaRqUbUFjgH+oEDRBPfJIUlFnh2K/1OkEBI+FCbBoECxQVQspqcnYSJvMkmRS0g0qOYr6aopzLz8wjHAP9yLxyMcQMiOI+tQP3t1ICuEUz2Sj+m0dNhdix0oSSNe21ZXTSETsUUMBYoCwjTuVZyw489vvAoJ835MdvXDLScJ93UixB9icgEMbEcq/bvKhoAC9v/BoUBxwa9qZuUkVqjkURF1p372atudh8OyshCSzwhxotqvl7b3WS7v0JoSDRQoinBphwhkC4eqOAiS4M2t7gWxqFBIEzIQFYu52Zoif4btSg90knVAVEAvDlItO5excuY5XsWAX3HipR55FdCnai9hPSXkHELk+0nWJsYHTmL1QwuKJtjZpxOz9UJ3KvyWncsyFhv5/yrP4nQeIWlDXt4x49Ru7JZ5ZLjkEw4UKA60rb0fAFC3qAkARQjpR3a0k98LE1H3/IoUIUror0KId6xEStva+zNjg2Do/qM5uV9XUqFAUUAIFUIEslAI2yJhFsZuszUnS0oS1shFp852RXIJ0eaEUKlb1MQ6HDIUKIT4JK5BXsWUbBf9Q0iaMVtBnMS8k6OsyrIPCQ6dZAnJMcJ0no0C89IpIUnEqu2YAydYh8OFAoWQHERFpNiJE/qhkLThlvzQCqt2QkfYaKFAISRPsdtDJAlOfLSikCiR20Fpe19W7hI7nJZ3SDRQoBBCCMlb7ERwEEuiECl0kg0XzwLllVdewZw5c1BZWYmCggI8//zzWccNw8CKFStQWVmJoUOHor6+HgcOHMg6p7e3FwsWLMCoUaMwbNgw3HzzzThy5EigL0JI2vDiiyISCCYxVL5+9mrUz15NawoJBSEivCbedBMwtKSEj2eB8umnn6Kurg5PPPGE5fFHH30Ua9aswRNPPIHXX38d5eXlmDlzJk6cOJE5p7GxEdu2bcOWLVvw6quv4uTJk7jppptw5swZ/9+EkBRiJTySKkSskDt5dvgkDIIIXyeRkittLJcpMAzD8P3hggJs27YNt9xyC4B+60llZSUaGxuxbFn/j9fb24uysjKsXr0a99xzD7q7u/G5z30Omzdvxu233w4A+Oijj1BVVYXt27fjxhtvdL1vT08PSkpK0N3djREjRvh9fELynvrZqxPbkVo5Lib1WUluYxYppe19npZ4kpA/KF/wMn5r9UE5ePAgOjo6MGvWrMx7RUVFuPbaa7F3714AQGtrKz777LOscyorK1FbW5s5hxCihyR3qGZTe5KfleQPfnenZ/2MHq0CpaOjAwBQVlaW9X5ZWVnmWEdHBwoLC1FaWmp7jpne3l709PRkvQghuY3sYMjOn4SF2XridXO/JES9pZVQongKCgqy/jYMY8B7ZpzOWbVqFUpKSjKvqqoqbc9KCCGEkOShVaCUl5cDwABLSGdnZ8aqUl5ejr6+PnR1ddmeY2b58uXo7u7OvA4fPqzzsQnJO+oWNWX5ePhJVBUGImJH/J+QqKEzdu6gVaBUV1ejvLwczc3Nmff6+vqwZ88eTJ8+HQAwefJkDBkyJOucY8eOYf/+/ZlzzBQVFWHEiBFZL0KIM101hVmCIG7MgkneoZmQOLBb5tlxZF3mReLD82aBJ0+exG9/+9vM3wcPHsTbb7+NkSNHYuzYsWhsbMTKlSsxYcIETJgwAStXrsRFF12Eb37zmwCAkpIS3H333ViyZAkuvvhijBw5EkuXLsXEiRMxY8YMfd+MEOK683FU2IkkihMSNm5Zi502DCTx4lmgvPHGG7juuusyfy9evBgAcMcdd2Djxo144IEHcOrUKXzve99DV1cXrr76auzatQvFxcWZzzQ1NWHw4MG47bbbcOrUKdxwww3YuHEjBg0apOErEULa1t6f6MRnFCaEEDcC5UGJC+ZBIcQdIVBK2/tiFQR1i5oGWHEoUEiUyGLdzqIoW1G4tBMeseVBIYQkB2HajluczLzrNRz66nljLcUJIUQFChRC8pg4NzMTs9Zdh2pQ/D67GpIb0HqSHDz7oBBCiBvy8hLahwGId5mJEJJ7cFpDCNFKljghJEa8OorTepIsaEEhhGgjyZFDJJ241UnhHEtxkjxoQSGEhA6Xd0jUCGFizOjCpiVrLM+hOEk2tKAQQrQhnHLlxGwUJyRpMDFbbsA8KIQQQvKSukVN2LRkDea1fRsFu0sBDPSNohUlWpgHhRBCCPFAw5iFcT8CMUGBQgghJC9pW3s/NndNx8mDJa7n2m0cSOKDPiiEEELyCjlyZ9OSvdiKqzJ/mzfQpDBJLrSgEEIIyUtOjD8LABhe3e16Lp25kwcFCiGEkNRglUCQ4iSZUKAQQgjJK5z2oOqqKURXTWGET0P8QoFCCCEk71DZKJOWk2RDgUIIISQvEbtozxrbnvW+vMxDkZJcKFAIIYTkNfNK9wLoFybcxDJ3YJgxIYSQvKRt7f1Y2vYamn8xFcDAEGMvyNs3ALS8RAEFCiGEkLykPx/KVO3XpTiJBgoUQggheUlpe19WxE6Q5R2KkuihDwohhJC8xCwqzOHF5mUbkiwoUAghhKQKWbhQpCQXChRCCCF5i5wPpW3t/fhk/qe4dN1j+GT+p0zalnAKDMMw4n4Ir/T09KCkpATd3d0YMWJE3I9DCCEkR2g7VIXNXdMBIBPdU9reRx+TiPAyftNJlhBCSGq49fnGzP+Lz/1LK0oyoUAhhBCSGj5YuORc+HE/KinxSTxQoBBCCEkVFCW5AZ1kCSGEEJI4KFAIIYQQkjgoUAghhBCSOChQCCGEEJI4KFAIIYQQkjgoUAghhBCSOChQCCGEEJI4KFAIIYQQkjgoUAghhBCSOChQCCGEEJI4KFAIIYQQkjgoUAghhBCSOChQCCGEEJI4KFAIIYQQkji0C5QVK1agoKAg61VeXp45bhgGVqxYgcrKSgwdOhT19fU4cOCA7scghBBCQqHtUFXmRcJjcBgXveKKK7B79+7M34MGDcr8/9FHH8WaNWuwceNGXHbZZfjRj36EmTNnor29HcXFxWE8DiGEEBIIipHoCUWgDB48OMtqIjAMA48//jgeeughzJ07FwDw1FNPoaysDM888wzuueeeMB6HEEII8YQQJHVjDw8QJ3VjD8fxSKkjFB+U9957D5WVlaiursY3vvENfPDBBwCAgwcPoqOjA7NmzcqcW1RUhGuvvRZ79+61vV5vby96enqyXoQQQkjYUJzEh3aBcvXVV2PTpk146aWX8POf/xwdHR2YPn06Pv74Y3R0dAAAysrKsj5TVlaWOWbFqlWrUFJSknlVVdHURrJxWhOW3+PaMSHEDfYPyUD7Ek9DQ0Pm/xMnTsS0adMwfvx4PPXUU5g6dSoAoKCgIOszhmEMeE9m+fLlWLx4cebvnp4eihQF5EaWJtXvJlLk91TKRfU83cgmZkIISRuh+KDIDBs2DBMnTsR7772HW265BQDQ0dGBioqKzDmdnZ0DrCoyRUVFKCoqCvtR84Ygg3GU2A3AXoVVWLMds+VF9Xn83sfu2kn87QjJR9z6ErbDaAldoPT29uI///M/8ed//ueorq5GeXk5mpub8aUvfQkA0NfXhz179mD16tVhP0oqcGpgSZ2RB3nmoOLEbn3Z7rqq69F+OjqVcnC6JyGE5BPaBcrSpUsxZ84cjB07Fp2dnfjRj36Enp4e3HHHHSgoKEBjYyNWrlyJCRMmYMKECVi5ciUuuugifPOb39T9KKkhDeulUX1Hr/fx+1xBvo/Xz1LQEDKQNPSbuY52gXLkyBH81V/9Ff7whz/gc5/7HKZOnYrXXnsN48aNAwA88MADOHXqFL73ve+hq6sLV199NXbt2sUcKB7QbTWQMYfUhTW4sXMghBDiRIFhGEbcD+GVnp4elJSUoLu7GyNGjIj7cSIhiQO6H/GSxO+RJmhNIWQgXvoltqFgeBm/Q/dBIf5J+mDuJQqGJAOV34IdMCEkCVCgRIguB8skYfWsHOByG0YNkXzEq0M8iR8KFISTKVCl0ufr4M4Gn/vka90k6UQ1Ks8NtoHg7Dv8ReVzc9oH5dX9lRhe3J8M12vYZhDywfJBiF/YSZNcQWefzHofDPFbnDxxFtfUfpQuH5QoxQGFCEkzzMlC0giXPtXRNUbmjUAhhESPW8g6IXHgtBMx0UMU5RrKbsaEEMJNGUncsP6FQ1TlSgsKIYQQ36gOVlFa1Gg50U8c5UmBQgghxBdeBi07H46gEWN2PlHy/3UNrmnzQ4lb5FGgEEIIiQTVAc8sOvx8jngjiWVHgUII0Yo8oKRpthk1dgOKzjJ3+x2jSHaWpIEzX+tzkspYhgKFEKIN0YHna0eeFJwGFL+DjdNvlralDSvy9fsnVZwAFCiEEE3kawceF1EPHHa+IOL9tOe/yQeRlmQxYgUFCiEkELneacdBrgwUulLE5wv5IFJyCQoUQgjRSFoH77SQi5akXK2TFCiEEF/kSuccBbk6AJBgcFPNcKFAIYQQD1CMECeSEMGWL3WUAoUQQs7BGTHRhbkuhV2P8kWUyFCgEEJ8YZ4pBu0gkyoE8rHjJ9Gj23clDfWSAoUQEgidacStCFO4pKGTJ8mD9U4NChRCSKJhZ05IOrkg7gcghBBCCDFDgUIIIYSQxEGBQgghhJDEQYFCCCGEkMRBgUIIIYSQxEGBQgghhJDEQYFCCCGEkMRBgUIIIYSQxEGBQgghhJDEQYFCCCGEkMRBgUIIIYSQxEGBQgghhJDEQYFCCCGEkMRBgUIIIYSQxEGBQgghhJDEQYFCCCGEkMRBgUIIIYSQxEGBQgghhJDEQYFCCCGEkMQRq0D56U9/iurqalx44YWYPHky/uM//iPOxyGEEEJIQohNoDz77LNobGzEQw89hLfeegt//ud/joaGBhw6dCiuRyKEEEJIQohNoKxZswZ33303/vZv/xZf/OIX8fjjj6Oqqgrr16+P65EIIYQQkhAGx3HTvr4+tLa24u/+7u+y3p81axb27t074Pze3l709vZm/u7u7gYAfHrybLgPSgghhBBtiHHbMAzXc2MRKH/4wx9w5swZlJWVZb1fVlaGjo6OAeevWrUKDz/88ID3b5w68FxCCCGEJJsTJ06gpKTE8ZxYBIqgoKAg62/DMAa8BwDLly/H4sWLM38fP34c48aNw6FDh1y/YFro6elBVVUVDh8+jBEjRsT9OImAZTIQlslAWCYDYZkMhGUyED9lYhgGTpw4gcrKStdzYxEoo0aNwqBBgwZYSzo7OwdYVQCgqKgIRUVFA94vKSlhRTExYsQIlokJlslAWCYDYZkMhGUyEJbJQLyWiaphIRYn2cLCQkyePBnNzc1Z7zc3N2P69OlxPBIhhBBCEkRsSzyLFy/GvHnzMGXKFEybNg0/+9nPcOjQIdx7771xPRIhhBBCEkJsAuX222/Hxx9/jB/+8Ic4duwYamtrsX37dowbN871s0VFRfjBD35gueyTVlgmA2GZDIRlMhCWyUBYJgNhmQwk7DIpMFRifQghhBBCIoR78RBCCCEkcVCgEEIIISRxUKAQQgghJHFQoBBCCCEkceSkQPnpT3+K6upqXHjhhZg8eTL+4z/+I+5HCo1XXnkFc+bMQWVlJQoKCvD8889nHTcMAytWrEBlZSWGDh2K+vp6HDhwIOuc3t5eLFiwAKNGjcKwYcNw880348iRIxF+C32sWrUKV155JYqLizF69GjccsstaG9vzzonbWWyfv16TJo0KZMsadq0adixY0fmeNrKw4pVq1ahoKAAjY2NmffSVi4rVqxAQUFB1qu8vDxzPG3lITh69Cj+5m/+BhdffDEuuugi/Omf/ilaW1szx9NWLn/yJ38yoJ4UFBRg/vz5ACIuDyPH2LJlizFkyBDj5z//ufHuu+8aixYtMoYNG2b87ne/i/vRQmH79u3GQw89ZDz33HMGAGPbtm1Zxx955BGjuLjYeO6554x9+/YZt99+u1FRUWH09PRkzrn33nuNSy65xGhubjbefPNN47rrrjPq6uqM06dPR/xtgnPjjTcaGzZsMPbv32+8/fbbxte+9jVj7NixxsmTJzPnpK1MXnjhBeNXv/qV0d7ebrS3txsPPvigMWTIEGP//v2GYaSvPMz85je/Mf7kT/7EmDRpkrFo0aLM+2krlx/84AfGFVdcYRw7dizz6uzszBxPW3kYhmF88sknxrhx44w777zT+H//7/8ZBw8eNHbv3m389re/zZyTtnLp7OzMqiPNzc0GAOPXv/61YRjRlkfOCZSrrrrKuPfee7Pe+8IXvmD83d/9XUxPFB1mgXL27FmjvLzceOSRRzLv/fGPfzRKSkqM//k//6dhGIZx/PhxY8iQIcaWLVsy5xw9etS44IILjJ07d0b27GHR2dlpADD27NljGAbLRFBaWmr8r//1v1JfHidOnDAmTJhgNDc3G9dee21GoKSxXH7wgx8YdXV1lsfSWB6GYRjLli0zrrnmGtvjaS0XmUWLFhnjx483zp49G3l55NQST19fH1pbWzFr1qys92fNmoW9e/fG9FTxcfDgQXR0dGSVR1FREa699tpMebS2tuKzzz7LOqeyshK1tbV5UWbd3d0AgJEjRwJgmZw5cwZbtmzBp59+imnTpqW+PObPn4+vfe1rmDFjRtb7aS2X9957D5WVlaiursY3vvENfPDBBwDSWx4vvPACpkyZgr/8y7/E6NGj8aUvfQk///nPM8fTWi6Cvr4+PP3007jrrrtQUFAQeXnklED5wx/+gDNnzgzYULCsrGzAxoNpQHxnp/Lo6OhAYWEhSktLbc/JVQzDwOLFi3HNNdegtrYWQHrLZN++fRg+fDiKiopw7733Ytu2bbj88stTWx4AsGXLFrz55ptYtWrVgGNpLJerr74amzZtwksvvYSf//zn6OjowPTp0/Hxxx+nsjwA4IMPPsD69esxYcIEvPTSS7j33nuxcOFCbNq0CUA664nM888/j+PHj+POO+8EEH15xJbqPggFBQVZfxuGMeC9NOGnPPKhzO677z688847ePXVVwccS1uZ1NTU4O2338bx48fx3HPP4Y477sCePXsyx9NWHocPH8aiRYuwa9cuXHjhhbbnpalcGhoaMv+fOHEipk2bhvHjx+Opp57C1KlTAaSrPADg7NmzmDJlClauXAkA+NKXvoQDBw5g/fr1+Na3vpU5L23lInjyySfR0NCAysrKrPejKo+csqCMGjUKgwYNGqDCOjs7Byi6NCA88J3Ko7y8HH19fejq6rI9JxdZsGABXnjhBfz617/GmDFjMu+ntUwKCwvx+c9/HlOmTMGqVatQV1eHtWvXprY8Wltb0dnZicmTJ2Pw4MEYPHgw9uzZg3Xr1mHw4MGZ75W2cpEZNmwYJk6ciPfeey+19aSiogKXX3551ntf/OIXcejQIQDp7U8A4He/+x12796Nv/3bv828F3V55JRAKSwsxOTJk9Hc3Jz1fnNzM6ZPnx7TU8VHdXU1ysvLs8qjr68Pe/bsyZTH5MmTMWTIkKxzjh07hv379+dkmRmGgfvuuw9bt27Fyy+/jOrq6qzjaSwTKwzDQG9vb2rL44YbbsC+ffvw9ttvZ15TpkzBX//1X+Ptt9/GpZdemspykent7cV//ud/oqKiIrX15Mtf/vKANAX/9V//ldm0Nq3lAgAbNmzA6NGj8bWvfS3zXuTl4cerN05EmPGTTz5pvPvuu0ZjY6MxbNgw48MPP4z70ULhxIkTxltvvWW89dZbBgBjzZo1xltvvZUJq37kkUeMkpISY+vWrca+ffuMv/qrv7IM+RozZoyxe/du48033zSuv/76nA2B++53v2uUlJQYLS0tWaFw//3f/505J21lsnz5cuOVV14xDh48aLzzzjvGgw8+aFxwwQXGrl27DMNIX3nYIUfxGEb6ymXJkiVGS0uL8cEHHxivvfaacdNNNxnFxcWZvjNt5WEY/SHogwcPNv7+7//eeO+994xf/vKXxkUXXWQ8/fTTmXPSWC5nzpwxxo4dayxbtmzAsSjLI+cEimEYxk9+8hNj3LhxRmFhofFnf/ZnmRDTfOTXv/61AWDA64477jAMoz8M7gc/+IFRXl5uFBUVGV/5yleMffv2ZV3j1KlTxn333WeMHDnSGDp0qHHTTTcZhw4diuHbBMeqLAAYGzZsyJyTtjK56667Mu3hc5/7nHHDDTdkxIlhpK887DALlLSVi8hXMWTIEKOystKYO3euceDAgczxtJWH4MUXXzRqa2uNoqIi4wtf+ILxs5/9LOt4GsvlpZdeMgAY7e3tA45FWR4FhmEYnm0/hBBCCCEhklM+KIQQQghJBxQohBBCCEkcFCiEEEIISRwUKIQQQghJHBQohBBCCEkcFCiEEEIISRwUKIQQQghJHBQohBBCCEkcFCiEEEIISRwUKIQQQghJHBQohBBCCEkcFCiEEEIISRz/PxbMLoEQjFtlAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "soilmvar = gfs.variables['Volumetric_Soil_Moisture_Content_depth_below_surface_layer']\n",
    "# flip the data in latitude so North Hemisphere is up on the plot\n",
    "soilm = soilmvar[0,0,::-1,:] \n",
    "print('shape=%s, type=%s, missing_value=%s' % \\\n",
    "      (soilm.shape, type(soilm), soilmvar.missing_value))\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "cs = plt.contourf(soilm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 32,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "##Packed integer data\n",
    "There is a similar feature for variables with `scale_factor` and `add_offset` attributes.\n",
    "\n",
    "- short integer data will automatically be returned as float data, with the scale and offset applied.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 32,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Dealing with dates and times\n",
    "- time variables usually measure relative to a fixed date using a certain calendar, with units specified like ***`hours since YY:MM:DD hh-mm-ss`***.\n",
    "- **`num2date`** and **`date2num`** convenience functions provided to convert between these numeric time coordinates and handy python datetime instances.  \n",
    "- **`date2index`** finds the time index corresponding to a datetime instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 34
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "name of time dimension = time1\n",
      "units = Hour since 2023-05-25T12:00:00Z, values = [  0.   3.   6.   9.  12.  15.  18.  21.  24.  27.  30.  33.  36.  39.\n",
      "  42.  45.  48.  51.  54.  57.  60.  63.  66.  69.  72.  75.  78.  81.\n",
      "  84.  87.  90.  93.  96.  99. 102. 105. 108. 111. 114. 117. 120. 123.\n",
      " 126. 129. 132. 135. 138. 141. 144. 147. 150. 153. 156. 159. 162. 165.\n",
      " 168. 171. 174. 177. 180. 183. 186. 189. 192. 195. 198. 201. 204. 207.\n",
      " 210. 213. 216. 219. 222. 225. 228. 231. 234. 237. 240. 243. 246. 249.\n",
      " 252. 255. 258. 261. 264. 267. 270. 273. 276. 279. 282. 285. 288. 291.\n",
      " 294. 297. 300. 303. 306. 309. 312. 315. 318. 321. 324. 327. 330. 333.\n",
      " 336. 339. 342. 345. 348. 351. 354. 357. 360. 363. 366. 369. 372. 375.\n",
      " 378. 381. 384.]\n"
     ]
    }
   ],
   "source": [
    "from netCDF4 import num2date, date2num, date2index\n",
    "timedim = sfctmp.dimensions[0] # time dim name\n",
    "print('name of time dimension = %s' % timedim)\n",
    "times = gfs.variables[timedim] # time coord var\n",
    "print('units = %s, values = %s' % (times.units, times[:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 35,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['2023-05-25 12:00:00', '2023-05-25 15:00:00', '2023-05-25 18:00:00', '2023-05-25 21:00:00', '2023-05-26 00:00:00', '2023-05-26 03:00:00', '2023-05-26 06:00:00', '2023-05-26 09:00:00', '2023-05-26 12:00:00', '2023-05-26 15:00:00']\n"
     ]
    }
   ],
   "source": [
    "dates = num2date(times[:], times.units)\n",
    "print([date.strftime('%Y-%m-%d %H:%M:%S') for date in dates[:10]]) # print only first ten..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 35,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "###Get index associated with a specified date, extract forecast data for that date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 37
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2023-05-28 15:57:27.760935\n",
      "index = 25, date = 2023-05-28 15:00:00\n"
     ]
    }
   ],
   "source": [
    "from datetime import datetime, timedelta\n",
    "date = datetime.now() + timedelta(days=3)\n",
    "print(date)\n",
    "ntime = date2index(date,times,select='nearest')\n",
    "print('index = %s, date = %s' % (ntime, dates[ntime]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 38
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "###Get temp forecast for Boulder (near 40N, -105W)\n",
    "- use function **`getcloses_ij`** we created before..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 39,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Boulder forecast valid at 2023-05-28 15:00:00 UTC = 297.6 K\n"
     ]
    }
   ],
   "source": [
    "lats, lons = gfs.variables['lat'][:], gfs.variables['lon'][:]\n",
    "# lats, lons are 1-d. Make them 2-d using numpy.meshgrid.\n",
    "lons, lats = np.meshgrid(lons,lats)\n",
    "j, i = getclosest_ij(lats,lons,40,-105)\n",
    "fcst_temp = sfctmp[ntime,j,i]\n",
    "print('Boulder forecast valid at %s UTC = %5.1f %s' % \\\n",
    "      (dates[ntime],fcst_temp,sfctmp.units))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 39,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##Simple multi-file aggregation\n",
    "\n",
    "What if you have a bunch of netcdf files, each with data for a different year, and you want to access all the data as if it were in one file?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 41
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-rw-r-- 1 8985332 May 17 15:27 data/prmsl.2000.nc\r\n",
      "-rw-rw-r-- 1 8968789 May 17 15:27 data/prmsl.2001.nc\r\n",
      "-rw-rw-r-- 1 8972796 May 17 15:27 data/prmsl.2002.nc\r\n",
      "-rw-rw-r-- 1 8974435 May 17 15:27 data/prmsl.2003.nc\r\n",
      "-rw-rw-r-- 1 8997438 May 17 15:27 data/prmsl.2004.nc\r\n",
      "-rw-rw-r-- 1 8976678 May 17 15:27 data/prmsl.2005.nc\r\n",
      "-rw-rw-r-- 1 8969714 May 17 15:27 data/prmsl.2006.nc\r\n",
      "-rw-rw-r-- 1 8974360 May 17 15:27 data/prmsl.2007.nc\r\n",
      "-rw-rw-r-- 1 8994260 May 17 15:27 data/prmsl.2008.nc\r\n",
      "-rw-rw-r-- 1 8974678 May 17 15:27 data/prmsl.2009.nc\r\n",
      "-rw-rw-r-- 1 8970732 May 17 15:27 data/prmsl.2010.nc\r\n",
      "-rw-rw-r-- 1 8976285 May 17 15:27 data/prmsl.2011.nc\r\n"
     ]
    }
   ],
   "source": [
    "!ls -ldgG data/prmsl*nc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 42
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**`MFDataset`** uses file globbing to patch together all the files into one big Dataset.\n",
    "You can also pass it a list of specific files.\n",
    "\n",
    "Limitations:\n",
    "\n",
    "- It can only  aggregate the data along the leftmost dimension of each variable.\n",
    "- only works with `NETCDF3`, or `NETCDF4_CLASSIC` formatted files.\n",
    "- kind of slow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 43,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting date = 2000-01-01 00:00:00\n",
      "ending date = 2011-12-31 00:00:00\n",
      "times shape = 4383\n",
      "prmsl dimensions = ('time', 'lat', 'lon'), prmsl shape = (4383, 91, 180)\n"
     ]
    }
   ],
   "source": [
    "mf = netCDF4.MFDataset('data/prmsl*nc')\n",
    "times = mf.variables['time']\n",
    "dates = num2date(times[:],times.units)\n",
    "print('starting date = %s' % dates[0])\n",
    "print('ending date = %s'% dates[-1])\n",
    "prmsl = mf.variables['prmsl']\n",
    "print('times shape = %s' % times.shape)\n",
    "print('prmsl dimensions = %s, prmsl shape = %s' %\\\n",
    "     (prmsl.dimensions, prmsl.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 43,
     "slide_type": "subslide"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Closing your netCDF file\n",
    "\n",
    "It's good to close netCDF files, but not actually necessary when Dataset is open for read access only.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 45
    },
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "f.close()\n",
    "gfs.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "internals": {
     "frag_helper": "fragment_end",
     "frag_number": 45,
     "slide_helper": "subslide_end"
    },
    "slide_helper": "slide_end",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "##That's it!\n",
    "\n",
    "Now you're ready to start exploring your data interactively.\n",
    "\n",
    "To be continued with **Writing netCDF data** ...."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
